Evaluating Short-Term Restoration Approaches for Degraded Grasslands in Kenya

Author

Ridhin Paul

Summary

This project evaluates the short-term effectiveness of restoration methods for degraded grasslands in Kenya. The treatments include Tilage (T), the addition of native soil inoculum (S) and Manure (M). “A grazing exclosure was installed at the campus of the Kabianga University in Kenya in April 2022, when restoration treatments were laid out and baseline soil samples taken and analysed for key soil physical, chemical and biological properties. In November 2022, fresh soil samples were taken to track short- term changes in the soil. In addition, we assessed the plant community species composition, and measured in situ soil respiration rates.” Refer “Grassland restoration experiment Kabianga Kenya.pdf” for more info.

Abrevations used

Variables and abbreviations Explanation
M Manure
S Soil inoculum
T Tillage
C Control
YIELD Yield in dry weight (kg/625 m^2)
DM Dry matter (%)
Tsoil Temperature of the topsoil (0C)
Tair Temperature of the air (0C)
MOIS Moisture of the topsoil (%)
pH pH in water (dimensionless)
EC Electrical Conductivity (µs m-1)
MBC Microbial Biomass Carbon (mg C kg soil-1)
MBN Microbial Biomass Nitrogen (mg N kg soil-1)
Nconc Soil total Nitrogen (%)
Cconc Soil total Carbon (%)
P Olsen P (μg P g-1 soil)
BD Bulk density (g/cm3)
NH4 Ammonium-Nitrogen (µg NH4-N g-1 DW)
NO3 Nitrate-Nitrogen (µg NO3-N g-1 DW)

Loading required packages

library(tidyverse)
library(readxl)
library(GGally)
library(VIM)
library(FactoMineR)
library(stats)
library(factoextra)
library(ggrepel)

Loading all datasets

Some unwanted columns have also been removed and -1 is replaced as 0 for easy interpretability.

#Kabianga
kab.nov.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                         sheet = 'Kab_nov_2022', 
                         na = c("NA", NA)) |>
  select(-c(...21,...22)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))
kab.dec.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                         sheet = 'Kab_dec_2023', 
                         na = c("NA", NA)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1)) 

#Thurgem
thur.jun.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx" , 
                          sheet = 'Thur_june_2022', 
                          na = c("NA", NA)) |>
  select(-c(...13,...14)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

thur.july.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                           sheet ='Thur_july_2023', 
                           na = c("NA", NA)) |>
  select(-c(...15,...16)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

thur.dec.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx" ,
                          sheet = 'Thur_dec_2023',
                          na = c("NA", NA)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

#Kapsarok
kaps.jun.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                          sheet = 'Kaps_june_2022', 
                          na = c("NA", NA)) |>
  select(-c(...13,...14))  |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

Data exploration

Missing value check

aggr(kab.nov.22, prop = F, numbers = T) 

aggr(kab.dec.23, prop = F, numbers = T)

aggr(thur.jun.22, prop = F, numbers = T)

aggr(thur.july.23, prop = F, numbers = T)

aggr(thur.dec.23, prop = F, numbers = T)

aggr(kaps.jun.22, prop = F, numbers = T)

For the dataset kab.nov.23 the plot indicate the last row has all NA entries which is removed. And for the dataset kaps.jun.22 there are 2 NA entries in the YIELD column, here the average based on the type of treatment for which the value is missing is taken and NA was replaced. Since there are only 2 NA values, code is not written for the entire column.

kab.nov.22 <- kab.nov.22|>
  na.omit()

mean_yield_1 <- kaps.jun.22 |>
  filter(Code == "S") |>
  pull(YIELD) |>
  mean(na.rm = TRUE)

mean_yield_2 <- kaps.jun.22 |>
  filter(Code == "ST") |>
  pull(YIELD) |>
  mean(na.rm = TRUE)

kaps.jun.22 <- kaps.jun.22 |>
  mutate(YIELD = if_else(Code == "S" & is.na(YIELD), mean_yield_1, YIELD),
         YIELD = if_else(Code == "ST" & is.na(YIELD), mean_yield_2, YIELD)) 

Checking data types

For the variables Plot No., YIELD, DM, SPE, Tsoil, Tair, MOIS, RES, pH, EC, MBC, MBN, Nconc, Cconc, P, BD, NH4, NO3, Block and for each treatment type M, S and T numeric values are expected as entry.

glimpse(kab.nov.22)
Rows: 32
Columns: 20
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Manure+Soil+Tillage", "Control", "Manure", "Manure…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ Tsoil      <dbl> 18.9, 18.6, 18.7, 18.7, 18.7, 18.7, 18.8, 19.7, 20.7, 21.8,…
$ Tair       <dbl> 19.7, 19.8, 20.2, 20.3, 20.6, 20.9, 21.7, 24.6, 29.9, 31.2,…
$ MOIS       <dbl> 53.6, 37.0, 51.7, 51.6, 40.4, 52.4, 37.6, 36.6, 40.4, 38.7,…
$ YIELD      <dbl> 21.2, 13.8, 15.2, 17.8, 16.0, 16.4, 18.1, 17.1, 16.8, 15.2,…
$ NH4        <dbl> 0.7, 1.1, 2.0, 1.6, 1.2, 1.8, 0.5, 0.7, 2.1, 0.5, 1.1, 1.2,…
$ NO3        <dbl> 22.8, 31.0, 23.6, 5.3, 12.6, 24.6, 25.1, 25.3, 45.6, 14.6, …
$ P          <dbl> 1.4, 2.5, 0.6, 1.4, 1.7, 1.3, 1.0, 2.1, 1.4, 1.7, 2.0, 3.2,…
$ pH         <dbl> 5.0, 4.9, 4.7, 4.5, 4.6, 4.7, 4.7, 4.7, 4.6, 4.6, 4.4, 4.6,…
$ EC         <dbl> 41.9, 47.4, 35.1, 11.6, 19.7, 38.5, 35.2, 36.8, 59.9, 22.0,…
$ MBC        <dbl> 461.6, 512.7, 577.2, 574.4, 440.7, 514.2, 440.1, 422.5, 514…
$ MBN        <dbl> 88.0, 109.9, 108.7, 109.4, 85.9, 93.4, 81.8, 83.8, 98.5, 93…
$ Nconc      <dbl> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.4, 0.5, 0.5, 0.5, 0.4,…
$ Cconc      <dbl> 5.1, 5.4, 5.4, 5.0, 5.1, 5.3, 4.9, 4.8, 5.2, 5.3, 5.8, 4.6,…
glimpse(kab.dec.23)
Rows: 32
Columns: 10
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Manure+Soil+Tillage", "Control", "Manure", "Manure…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ BD         <dbl> 0.3589024, 0.4051220, 0.3829878, 0.4227439, 0.4143293, 0.38…
$ P          <dbl> 5.9156125, 3.1895944, 1.5363629, 1.3363361, 1.8315103, 1.40…
$ YIELD      <dbl> 10.43400, 6.74075, 6.88250, 6.60400, 5.60350, 6.98075, 11.3…
glimpse(thur.jun.22)
Rows: 32
Columns: 12
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Ti…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ Tsoil      <dbl> 21.8, 22.5, 22.8, 23.5, 24.0, 24.4, 24.8, 25.7, 26.6, 27.7,…
$ Tair       <dbl> 24.1, 25.7, 26.7, 27.4, 28.0, 27.8, 28.7, 30.7, 32.8, 33.8,…
$ MOIS       <dbl> 44.7, 46.3, 45.0, 42.7, 28.0, 43.2, 33.3, 38.1, 33.0, 39.5,…
$ YIELD      <dbl> 4.2, 5.8, 3.8, 3.1, 3.5, 5.1, 3.9, 2.5, 3.9, 4.2, 1.6, 6.7,…
$ EC         <dbl> 136.0, 149.5, 165.3, 182.6, 176.0, 190.4, 152.0, 148.0, 128…
glimpse(thur.july.23)
Rows: 32
Columns: 14
$ Plot_No   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ Block     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, …
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Til…
$ Code      <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "M…
$ M         <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, …
$ S         <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, …
$ T         <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, …
$ Tsoil     <dbl> 23.7, 23.6, 23.9, 24.8, 25.1, 25.0, 25.2, 26.7, 26.4, 28.0, …
$ Tair      <dbl> 25.3, 26.3, 28.0, 30.4, 32.3, 32.1, 32.4, 35.3, 33.5, 33.9, …
$ MOIS      <dbl> 48.1, 52.3, 47.3, 40.1, 47.9, 39.6, 47.8, 45.9, 51.6, 44.4, …
$ DW        <dbl> 162.7, 214.9, 270.3, 262.5, 298.1, 648.9, 646.1, 465.1, 253.…
$ YIELD     <dbl> 4.1, 5.4, 6.8, 6.6, 7.5, 16.2, 16.2, 11.6, 6.3, 8.4, 4.0, 19…
$ Nconc     <dbl> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, …
$ Cconc     <dbl> 1.0, 1.1, 1.3, 1.0, 1.0, 1.4, 1.2, 1.2, 0.9, 0.9, 0.8, 1.4, …
glimpse(thur.dec.23)
Rows: 32
Columns: 10
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Ti…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ MOIS       <dbl> 55.30603, 63.17742, 60.24427, 56.02128, 59.56154, 60.94444,…
$ P          <dbl> 2.7800613, 1.6521765, 1.6484200, 1.5918039, 1.6554478, 1.41…
$ YIELD      <dbl> 2.59225, 2.85375, 1.30200, 1.55025, 2.62850, 4.39375, 4.625…
glimpse(kaps.jun.22)
Rows: 32
Columns: 12
$ Plot_No.  <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ Block     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, …
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Til…
$ Code      <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "M…
$ M         <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, …
$ S         <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, …
$ T         <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, …
$ Tsoil     <dbl> 21.9, 21.9, 22.4, 23.5, 24.2, 24.7, 25.4, 26.2, 26.7, 28.7, …
$ Tair      <dbl> 24.3, 24.9, 26.1, 27.0, 28.8, 30.9, 32.5, 31.4, 32.0, 35.6, …
$ MOIS      <dbl> 12.0, 15.1, 19.9, 19.3, 25.4, 28.8, 26.2, 22.4, 18.9, 24.3, …
$ YIELD     <dbl> 1.866667, 2.900000, 1.400000, 2.300000, 10.300000, 6.500000,…
$ EC        <dbl> 106.0, 136.6, 140.8, 156.4, 110.8, 127.8, 101.8, 184.4, 140.…

All results were as expected.

Check for multicolinearity

  1. Kabianga November 2022

    kab.nov.22 |>
      select(-`Plot No.`, - Treatment, -M, -S, -T) |> 
      ggpairs()

    For the dataset a high correlation is observed between,

    Var 1 Var 2 r
    Tair Tsoil 0.909
    Nconc Tair 0.615
    YIELD MOIS 0.571
    MBN NH4 0.651
    EC NO3 0.980
    Nconc NO3 0.603
    MBN MBC 0.891
    Cconc MBC 0.627
    Cconc Nconc 0.731
  1. Kabianga December 2023

    kab.dec.23 |>
      select(-`Plot No.`, - Treatment, -M, -S, -T) |> 
      ggpairs()

    No significant correlations were seen for the dataset.

  1. Thurgem June 2022

    thur.jun.22 |>
      select(-`Plot No.`, - Treatment, -M, -S, -T) |> 
      ggpairs()

    For the dataset a high correlation is observed between,

    Var 1 Var 2 r
    Tair Tsoil 0.974
  1. Thurgem July 2023

    thur.july.23 |>
      select(-Plot_No, - Treatment, -M, -S, -T) |> 
      ggpairs()

    For the dataset a high correlation is observed between,

    Var 1 Var 2 r
    Tair Tsoil 0.912
    YIELD DW 1
  1. Thurgem Dec 2023

    thur.dec.23 |>
      select(- `Plot No.`, - Treatment, -M, -S, -T) |> 
      ggpairs()

    No significant correlations were observed between variables.

  1. Kapsorok June 2022

    kaps.jun.22 |>
      select(-Plot_No., - Treatment, -M, -S, -T) |> 
      ggpairs()

    For the dataset a high correlation is observed between,

    Var 1 Var 2 r
    Tair Tsoil 0.962

Correlation between Tair and Tsoil is consistent across all datasets.

Histogram of YIELD by Treatment faceted by Block

  1. Kabianga November 2022

    ggplot(kab.nov.22, aes(x = Code, y = YIELD)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Histogram of YIELD by Treaetment faceted by Block (Kabianga Nov 2022)", x
           = "Treatment", y = "Yield") +
        facet_wrap(~Block) +
      theme_minimal() 

  1. Kabianga December 2023

    ggplot(kab.dec.23, aes(x = Code, y = YIELD)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Histogram of YIELD by Treaetment faceted by Block (Kabianga Dec 2023)", x
           = "Treatment", y = "Yield") +
        facet_wrap(~Block) +
      theme_minimal() 

  1. Thurgem June 2022

    ggplot(thur.jun.22, aes(x = Code, y = YIELD)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem June 2022)", x
           = "Treatment", y = "Yield") +
        facet_wrap(~Block) +
      theme_minimal() 

  1. Thurgem July 2023

    ggplot(thur.july.23, aes(x = Code, y = YIELD)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem July 2023)", x
           = "Treatment", y = "Yield") +
        facet_wrap(~Block) +
      theme_minimal() 

  1. Thurgem Dec 2023

    ggplot(thur.dec.23, aes(x = Code, y = YIELD)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem Dec 2023)", x =
             "Treatment", y = "Yield") +
        facet_wrap(~Block) +
      theme_minimal() 

  1. Kapsorok June 2022

    ggplot(kaps.jun.22, aes(x = Code, y = YIELD)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of YIELD by Treaetment faceted by Block (Kapsorok June 2022)", x = "Treatment", y = "Yield") +
        facet_wrap(~Block) +
      theme_minimal() 

For all dataset the variation of YIELD across treatment across each block seem to vary slightly however it’s statistical significance need to be checked for each dataset.

Average yield across treatment

  1. Kabianga November 2022

    kab.nov.22 |>
      summarise(m_y_t = mean(YIELD), .by = Code) |>
      ggplot(aes(x = Code, y = m_y_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean YIELD by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean Yield") +
      theme_minimal()

For the dataset, it can be observed that for treatments MST, MT, ST and T there is a reduction in YIELD compared to C (Control).

  1. Kabianga December 2023

    kab.dec.23 |>
      summarise(m_y_t = mean(YIELD), .by = Code) |>
      ggplot(aes(x = Code, y = m_y_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean YIELD by Treaetment (Kabianga December 2023)", x = "Treatment", y = "Mean
           Yield") +
      theme_minimal()

For the dataset, it can be observed that for treatments MSTand M there is a reduction in YIELD compared to C (Control) and T the increase is fairly low.

  1. Thurgem June 2022

    thur.jun.22 |>
      summarise(m_y_t = mean(YIELD), .by = Code) |>
      ggplot(aes(x = Code, y = m_y_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem June 2022)", x = "Treatment", y = "Mean Yield") +
      theme_minimal()

For the dataset, it can be observed that for treatments M, MS, MST, S, MT, ST and T there is a reduction in YIELD compared to C (Control).

  1. Thurgem July 2023

    thur.july.23 |>
      summarise(m_y_t = mean(YIELD), .by = Code) |>
      ggplot(aes(x = Code, y = m_y_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem July 2023)", x = "Treatment", y = "Mean Yield") +
      theme_minimal()

For the dataset, it can be observed that for treatments M, MST, and ST there is a reduction in YIELD compared to C (Control) and for T a significant increase can be seen.

  1. Thurgem Dec 2023

    thur.dec.23 |>
      summarise(m_y_t = mean(YIELD), .by = Code) |>
      ggplot(aes(x = Code, y = m_y_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem Dec 2023)", x = "Treatment", y = "Mean Yield") +
      theme_minimal()

For the dataset, it can be observed that for every treatment there is an increase in YIELD compared to treatment.

  1. Kapsorok June 2022

    kaps.jun.22 |>
      summarise(m_y_t = mean(YIELD), .by = Code) |>
      ggplot(aes(x = Code, y = m_y_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean YIELD by Treaetment (Kapsorok June 2022)", x = "Treatment", y = "Mean Yield") +
      theme_minimal()

For the dataset, it can be observed that for treatments M, MST and S there is a reduction in YIELD compared to C (Control) and for the rest of the treatments there is a significant increase in YIELD.

Visually tillage seem to have a negative effect on YIELD and manure combined with soil inoculum seems to have a slight positive effect on productivity but cannot be concluded.

Average MBC and MBN across treatment

  1. Kabianga November 2023

    It seem that manure slightly increase both MBC and MBN which is expected because MBC and MBN had a high corelation. It is also observed that tillage has a slightly negative effect.

    kab.nov.22 |>
      summarise(m_mbc_t = mean(MBC), .by = Code) |>
      ggplot(aes(x = Code, y = m_mbc_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean MBC by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean MBC") +
      theme_minimal()

    kab.nov.22 |>
      summarise(m_mbn_t = mean(MBN), .by = Code) |>
      ggplot(aes(x = Code, y = m_mbn_t)) +
      geom_bar(color = "black", alpha = 0.7, stat = "identity") +
      labs(title = "Bar chart of mean MBN by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean MBN") +
      theme_minimal()

Primary analysis

Common function for brute force algorithm based on AIC. This was not implemented but only tested except for the dataset Kabianga November 2022 with dependent as YIELD because it yielded a sligtly better model. Hybrid stepwise algorithm is used based on AIC keeping the interaction effects as fixed for rest of the models. And for all models p<0.10 is considered significant.

#Pass the base model as the formulae and scope as a vector of variable names under consideration

bestModel <- function(base_model, scope, dataset)
  {
    all_combs <- list()
    #seq_along creates n sequence of elements from the inputed vector
    for (i in seq_along(scope)) 
      {
          #combn creates all possible n combinations, simplify is set to false to                 obtain the result as list
          combs <- combn(scope, i, simplify = FALSE)
          
          all_combs <- c(all_combs, combs)
      }
    
    combine <- function(all_combs) 
      {
        #deparse convert type formula to string
        paste(deparse(base_model), paste(all_combs, collapse =
                                   " + "), sep = " + ")
      }
    
    #applies the function to list all_combs and returns a list which is converted to vector
    formulas <- unlist(lapply(all_combs, combine))
    
    fitted_model <- list()
    
    for (formula in formulas)
      {
        model <- lm(as.formula(formula), data = dataset)
        
        #adjust below to consider other criteria apart from AIC and also ajust the following           logic accordingly after iteration
        fitted_model[[formula]] <- list(models = model, AIC = AIC(model))
  
      }
    
    #takes in list, applies a function and returns vector
    aic_values <- sapply(fitted_model, function(x) x$AIC)
    

    #returns the fitted model and corresponding AIC
    return(fitted_model[[which.min(aic_values)]])
  
  }
  1. Kabianga November 2022

    YIELD as the dependent variable.

    im <- lm(YIELD ~ T * S * M, data = kab.nov.22)
    
    model_kab_Nov22_YIELD <- step(im, 
                                  scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + Tsoil + Tair + MOIS + NH4 + NO3 +
                            P + pH + EC + MBC + MBN + Nconc + Cconc 
                        ),
                        direction = "both")  
    Start:  AIC=78.57
    YIELD ~ T * S * M
    
            Df Sum of Sq    RSS    AIC
    + MOIS   1   25.6598 200.42 76.710
    + MBN    1   15.9946 210.09 78.217
    <none>               226.08 78.565
    + NH4    1   11.2688 214.81 78.929
    + pH     1   10.2982 215.78 79.073
    + MBC    1    9.2920 216.79 79.222
    + Tsoil  1    7.6671 218.41 79.461
    + Tair   1    7.2481 218.83 79.523
    + EC     1    6.5224 219.56 79.628
    + NO3    1    5.6836 220.40 79.751
    + P      1    2.5974 223.49 80.195
    + Cconc  1    1.7285 224.35 80.320
    + Nconc  1    0.0001 226.08 80.565
    
    Step:  AIC=76.71
    YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq    RSS    AIC
    + pH     1   13.4101 187.01 76.494
    + MBN    1   13.0343 187.39 76.558
    <none>               200.42 76.710
    + P      1   10.0996 190.32 77.056
    + Tair   1    9.3024 191.12 77.189
    + NH4    1    9.1969 191.23 77.207
    + MBC    1    6.5987 193.82 77.639
    + Nconc  1    4.2286 196.19 78.028
    + Tsoil  1    3.4636 196.96 78.152
    - MOIS   1   25.6598 226.08 78.565
    + Cconc  1    0.0500 200.37 78.702
    + NO3    1    0.0310 200.39 78.705
    + EC     1    0.0071 200.42 78.709
    
    Step:  AIC=76.49
    YIELD ~ T + S + M + MOIS + pH + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq    RSS    AIC
    <none>               187.01 76.494
    + Tair   1   11.2434 175.77 76.510
    - pH     1   13.4101 200.42 76.710
    + Tsoil  1    8.3602 178.65 77.031
    + MBN    1    7.6768 179.34 77.153
    + NH4    1    5.8166 181.20 77.483
    + Nconc  1    4.1696 182.84 77.773
    + P      1    4.1162 182.90 77.782
    + MBC    1    2.3911 184.62 78.082
    + Cconc  1    0.0688 186.94 78.482
    + EC     1    0.0285 186.98 78.489
    + NO3    1    0.0031 187.01 78.494
    - MOIS   1   28.7718 215.78 79.073
    summary(model_kab_Nov22_YIELD)
    
    Call:
    lm(formula = YIELD ~ T + S + M + MOIS + pH + T:S + T:M + S:M + 
        T:S:M, data = kab.nov.22)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -6.9018 -1.4359  0.3215  1.4152  7.5303 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
    (Intercept)  18.6138    14.5300   1.281   0.2135  
    T             1.2088     2.8672   0.422   0.6774  
    S             5.0615     2.2751   2.225   0.0367 *
    M             1.7076     2.0662   0.826   0.4174  
    MOIS          0.2926     0.1590   1.840   0.0793 .
    pH           -3.4092     2.7143  -1.256   0.2223  
    T:S          -3.2100     3.0478  -1.053   0.3037  
    T:M          -2.1374     2.9200  -0.732   0.4719  
    S:M          -6.6686     3.0928  -2.156   0.0423 *
    T:S:M         5.2884     4.3188   1.225   0.2337  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 2.916 on 22 degrees of freedom
    Multiple R-squared:  0.4955,    Adjusted R-squared:  0.2891 
    F-statistic: 2.401 on 9 and 22 DF,  p-value: 0.04534
    #implementation of brute force algorithm
    scope = c("Tsoil", "Tair", "MOIS", "NH4", "NO3", "P", "pH", "EC", "MBC", "MBN", "Nconc", "Cconc")
    
    base_model = YIELD ~ T * S * M
    
    #run only if necessar of running bcz it takes a quite a bit of time the model m below is the result of this function run.
    # m <- bestModel(base_model, scope, kab.nov.22)
    
    # summary(m$model)
    
    #model m below is the result of this function run (best_Model)
    m <- lm(YIELD ~ T * S * M + MOIS + NH4 + MBN + Nconc, data = kab.nov.22) 
    summary(m)
    
    Call:
    lm(formula = YIELD ~ T * S * M + MOIS + NH4 + MBN + Nconc, data = kab.nov.22)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -4.3236 -1.3683 -0.1658  1.0606  4.9609 
    
    Coefficients:
                 Estimate Std. Error t value Pr(>|t|)   
    (Intercept)  -1.38995    8.65073  -0.161  0.87396   
    T             3.71723    2.81815   1.319  0.20206   
    S             6.29647    2.10707   2.988  0.00727 **
    M             1.82889    1.95230   0.937  0.36004   
    MOIS          0.42629    0.16314   2.613  0.01665 * 
    NH4           1.79159    1.49889   1.195  0.24596   
    MBN           0.12320    0.06815   1.808  0.08572 . 
    Nconc       -37.25304   15.70888  -2.371  0.02787 * 
    T:S          -2.37206    2.73333  -0.868  0.39578   
    T:M          -2.22589    2.72349  -0.817  0.42339   
    S:M          -7.41175    2.84954  -2.601  0.01709 * 
    T:S:M         3.27673    3.82547   0.857  0.40184   
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 2.695 on 20 degrees of freedom
    Multiple R-squared:  0.6081,    Adjusted R-squared:  0.3926 
    F-statistic: 2.821 on 11 and 20 DF,  p-value: 0.0212

    These results indicate that soil inoculum and the soil incolum combined with manure, moisture, MBN, and nitrogen concentration significantly influence yield, while other factors and interactions have less impact. Soil inoculum had a significant positive impact while soil incolum combined with manure had a negative impact.

    MBC as the dependent variable

    im <- lm(MBC ~ T * S * M, data = kab.nov.22)
    
    model_kab_Nov22_MBC <- step(im, 
                                  scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + Tsoil + Tair + MOIS + NH4 + NO3 +
                            P + pH + EC  + MBN + Nconc + Cconc 
                        ),
                        direction = "both")  
    Start:  AIC=261.84
    MBC ~ T * S * M
    
            Df Sum of Sq   RSS    AIC
    + MBN    1     55169 14274 213.21
    + Nconc  1     26477 42966 248.48
    + Cconc  1     23970 45473 250.29
    + NH4    1     16758 52685 255.00
    + EC     1     14522 54921 256.33
    + Tair   1     13595 55848 256.87
    + Tsoil  1     13523 55920 256.91
    + NO3    1     13448 55995 256.95
    + pH     1      5741 63702 261.08
    <none>               69443 261.84
    + P      1      3725 65718 262.08
    + MOIS   1       654 68789 263.54
    
    Step:  AIC=213.21
    MBC ~ T + S + M + MBN + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq   RSS    AIC
    + P      1      2644 11629 208.66
    + Cconc  1      1248 13026 212.29
    <none>               14274 213.21
    + NO3    1       819 13455 213.32
    + EC     1       533 13740 214.00
    + Nconc  1       519 13754 214.03
    + pH     1       349 13925 214.42
    + Tair   1       234 14040 214.69
    + NH4    1       119 14154 214.94
    + Tsoil  1        75 14199 215.05
    + MOIS   1        49 14224 215.10
    - MBN    1     55169 69443 261.84
    
    Step:  AIC=208.66
    MBC ~ T + S + M + MBN + P + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq   RSS    AIC
    + pH     1      1704  9925 205.59
    <none>               11629 208.66
    + NH4    1       475 11154 209.32
    + Cconc  1       312 11318 209.79
    + Tair   1        68 11561 210.47
    + MOIS   1        60 11569 210.49
    + Nconc  1        22 11607 210.60
    + Tsoil  1         7 11622 210.64
    + EC     1         4 11625 210.65
    + NO3    1         0 11629 210.66
    - P      1      2644 14274 213.21
    - MBN    1     54089 65718 262.08
    
    Step:  AIC=205.59
    MBC ~ T + S + M + MBN + P + pH + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq   RSS    AIC
    + NH4    1       741  9184 205.10
    <none>                9925 205.59
    + Cconc  1       382  9543 206.33
    + Tsoil  1       231  9694 206.83
    + Tair   1       219  9706 206.87
    + MOIS   1        58  9867 207.40
    + Nconc  1         5  9920 207.57
    + EC     1         3  9922 207.58
    + NO3    1         1  9924 207.58
    - pH     1      1704 11629 208.66
    - P      1      4000 13925 214.42
    - MBN    1     44873 54798 258.26
    
    Step:  AIC=205.1
    MBC ~ T + S + M + MBN + P + pH + NH4 + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq   RSS    AIC
    <none>                9184 205.10
    + Cconc  1       487  8697 205.36
    - NH4    1       741  9925 205.59
    + Tair   1       225  8959 206.31
    + Tsoil  1       160  9024 206.54
    + MOIS   1        74  9111 206.85
    + Nconc  1        26  9159 207.01
    + NO3    1        14  9170 207.06
    + EC     1        10  9174 207.07
    - pH     1      1970 11154 209.32
    - P      1      4606 13790 216.11
    - MBN    1     36389 45573 254.36
    summary(model_kab_Nov22_MBC)
    
    Call:
    lm(formula = MBC ~ T + S + M + MBN + P + pH + NH4 + T:S + T:M + 
        S:M + T:S:M, data = kab.nov.22)
    
    Residuals:
       Min     1Q Median     3Q    Max 
    -29.70 -10.73  -3.20  10.10  40.76 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept) 368.5011   126.1290   2.922  0.00844 ** 
    T           -21.6695    16.0188  -1.353  0.19123    
    S            10.4362    16.8553   0.619  0.54280    
    M            -4.2831    15.4943  -0.276  0.78505    
    MBN           4.3856     0.4927   8.902 2.15e-08 ***
    P           -22.8369     7.2111  -3.167  0.00485 ** 
    pH          -46.2661    22.3393  -2.071  0.05151 .  
    NH4         -14.6657    11.5453  -1.270  0.21856    
    T:S         -36.6072    24.2079  -1.512  0.14612    
    T:M           2.4122    21.8744   0.110  0.91329    
    S:M         -17.4707    23.1437  -0.755  0.45912    
    T:S:M        41.7015    33.6025   1.241  0.22896    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 21.43 on 20 degrees of freedom
    Multiple R-squared:  0.9063,    Adjusted R-squared:  0.8548 
    F-statistic: 17.59 on 11 and 20 DF,  p-value: 6.145e-08

    MBN as the dependent variable

    im <- lm(MBN ~ T * S * M, data = kab.nov.22)
    
    model_kab_Nov22_MBN <- step(im, 
                                  scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + Block + Tsoil + Tair + MOIS + NH4 + NO3 +
                            P + pH + EC  + MBC + Nconc + Cconc 
                        ),
                        direction = "both")  
    Start:  AIC=161.23
    MBN ~ T * S * M
    
            Df Sum of Sq     RSS    AIC
    + MBC    1   2378.00  615.24 112.60
    + Nconc  1   1137.18 1856.06 147.94
    + Block  1   1051.14 1942.10 149.38
    + NH4    1   1037.41 1955.83 149.61
    + Cconc  1    846.47 2146.77 152.59
    + Tsoil  1    640.22 2353.02 155.53
    + Tair   1    573.91 2419.34 156.42
    + EC     1    538.12 2455.12 156.89
    + NO3    1    434.67 2558.57 158.21
    <none>               2993.24 161.23
    + pH     1    180.42 2812.82 161.24
    + MOIS   1     18.69 2974.55 163.03
    + P      1      5.05 2988.19 163.17
    
    Step:  AIC=112.6
    MBN ~ T + S + M + MBC + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq     RSS    AIC
    + NH4    1     89.79  525.45 109.55
    + P      1     86.46  528.79 109.75
    <none>                615.24 112.60
    + Block  1     28.78  586.47 113.07
    + Nconc  1     21.08  594.17 113.49
    + Tsoil  1     17.77  597.47 113.66
    + Tair   1      7.04  608.20 114.23
    + EC     1      1.02  614.23 114.55
    + NO3    1      0.46  614.78 114.58
    + pH     1      0.38  614.87 114.58
    + Cconc  1      0.30  614.94 114.58
    + MOIS   1      0.17  615.08 114.59
    - MBC    1   2378.00 2993.24 161.23
    
    Step:  AIC=109.55
    MBN ~ T + S + M + MBC + NH4 + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq     RSS    AIC
    + P      1    103.39  422.07 104.54
    <none>                525.45 109.55
    + NO3    1     10.27  515.19 110.92
    + Block  1      8.73  516.72 111.02
    + Tsoil  1      4.64  520.81 111.27
    + EC     1      2.46  522.99 111.40
    + Tair   1      1.68  523.78 111.45
    + Cconc  1      1.40  524.05 111.47
    + Nconc  1      1.13  524.32 111.48
    + pH     1      0.97  524.49 111.49
    + MOIS   1      0.36  525.09 111.53
    - NH4    1     89.79  615.24 112.60
    - MBC    1   1430.38 1955.83 149.61
    
    Step:  AIC=104.54
    MBN ~ T + S + M + MBC + NH4 + P + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq     RSS    AIC
    + pH     1     40.79  381.28 103.29
    + Nconc  1     28.27  393.80 104.32
    + Tair   1     26.83  395.24 104.44
    <none>                422.07 104.54
    + Block  1     19.79  402.27 105.00
    + Tsoil  1     10.84  411.23 105.71
    + EC     1      7.14  414.92 106.00
    + MOIS   1      4.89  417.17 106.17
    + NO3    1      3.66  418.40 106.26
    + Cconc  1      1.65  420.41 106.42
    - P      1    103.39  525.45 109.55
    - NH4    1    106.72  528.79 109.75
    - MBC    1   1518.73 1940.80 151.36
    
    Step:  AIC=103.29
    MBN ~ T + S + M + MBC + NH4 + P + pH + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq     RSS    AIC
    + Tair   1     35.09  346.19 102.20
    + Tsoil  1     24.36  356.91 103.18
    + Block  1     23.93  357.35 103.22
    <none>                381.28 103.29
    + Nconc  1     17.90  363.38 103.75
    - pH     1     40.79  422.07 104.54
    + MOIS   1      4.98  376.30 104.87
    + EC     1      4.71  376.57 104.89
    + NO3    1      2.85  378.43 105.05
    + Cconc  1      0.06  381.21 105.28
    - NH4    1    117.87  499.14 109.91
    - P      1    143.21  524.49 111.49
    - MBC    1   1510.65 1891.93 152.55
    
    Step:  AIC=102.2
    MBN ~ T + S + M + MBC + NH4 + P + pH + Tair + T:S + T:M + S:M + 
        T:S:M
    
            Df Sum of Sq     RSS    AIC
    <none>                346.19 102.20
    - Tair   1     35.09  381.28 103.29
    + MOIS   1      7.22  338.97 103.53
    + Nconc  1      3.01  343.18 103.92
    + Tsoil  1      2.08  344.11 104.01
    + Block  1      1.84  344.35 104.03
    + NO3    1      1.74  344.45 104.04
    + EC     1      1.38  344.81 104.07
    + Cconc  1      0.21  345.98 104.18
    - pH     1     49.05  395.24 104.44
    - NH4    1    102.43  448.62 108.49
    - P      1    176.78  522.97 113.40
    - MBC    1   1297.49 1643.68 150.05
    summary(model_kab_Nov22_MBN)
    
    Call:
    lm(formula = MBN ~ T + S + M + MBC + NH4 + P + pH + Tair + T:S + 
        T:M + S:M + T:S:M, data = kab.nov.22)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -6.6871 -2.5162  0.5216  2.5615  5.6404 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept) -39.0872    28.6689  -1.363   0.1887    
    T             2.3904     3.3107   0.722   0.4791    
    S            -2.6106     3.4027  -0.767   0.4524    
    M             2.5024     3.0805   0.812   0.4267    
    MBC           0.1746     0.0207   8.439 7.52e-08 ***
    NH4           4.9858     2.1028   2.371   0.0285 *  
    P             5.0085     1.6079   3.115   0.0057 ** 
    pH            7.6942     4.6897   1.641   0.1173    
    Tair         -0.2306     0.1661  -1.388   0.1813    
    T:S           8.2018     4.9141   1.669   0.1115    
    T:M          -1.4562     4.3471  -0.335   0.7413    
    S:M           2.1094     4.6662   0.452   0.6563    
    T:S:M        -7.2801     6.8728  -1.059   0.3028    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 4.269 on 19 degrees of freedom
    Multiple R-squared:  0.909, Adjusted R-squared:  0.8515 
    F-statistic: 15.81 on 12 and 19 DF,  p-value: 2.073e-07

    MBC and MBN did not show any significant effects on treatment, but both had a very high significant effect on Phosphorus . MBC had a negative effect, whereas MBN had a positive effect on P.

  2. Kabianga December 2023

    im <- lm(YIELD ~ T * S * M, data = kab.dec.23)
    
    model_kab_Dec23_YIELD <- step(im, 
                                  scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + P + BD 
                        ),
                        direction = "both")  
    Start:  AIC=59.16
    YIELD ~ T * S * M
    
           Df Sum of Sq    RSS    AIC
    <none>              123.28 59.159
    + P     1   2.30009 120.98 60.556
    + BD    1   0.36336 122.92 61.064
    summary(model_kab_Dec23_YIELD)
    
    Call:
    lm(formula = YIELD ~ T * S * M, data = kab.dec.23)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -3.6739 -1.2265 -0.0384  1.1020  3.6000 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)   8.1330     1.1332   7.177 2.04e-07 ***
    T             0.6874     1.6026   0.429   0.6718    
    S             3.1477     1.6026   1.964   0.0612 .  
    M            -0.7659     1.6026  -0.478   0.6370    
    T:S          -1.2792     2.2664  -0.564   0.5777    
    T:M           1.2228     2.2664   0.540   0.5945    
    S:M          -1.2809     2.2664  -0.565   0.5772    
    T:S:M        -2.4089     3.2052  -0.752   0.4596    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 2.266 on 24 degrees of freedom
    Multiple R-squared:  0.3125,    Adjusted R-squared:  0.112 
    F-statistic: 1.559 on 7 and 24 DF,  p-value: 0.1958

    Soil inoculum showed a positive effect, which was significant and soil inoculum combined with manure had a negative effect, even though not statistically significant. The effect of treatment seems comparable to the initial measurement, even though a reduction in YIELD is observed.

  3. Thurgem June 2022

    im <- lm(YIELD ~ T * S * M, data = thur.jun.22)
    
    model_Thur_Jun22 <- step(im,
                            scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + Tair + Tsoil + MOIS + EC 
                        ),
                        direction = "both")  
    Start:  AIC=26.79
    YIELD ~ T * S * M
    
            Df Sum of Sq    RSS    AIC
    + MOIS   1   11.4629 33.370 19.341
    <none>               44.832 26.790
    + Tair   1    0.2215 44.611 28.632
    + EC     1    0.1946 44.638 28.651
    + Tsoil  1    0.1631 44.669 28.674
    
    Step:  AIC=19.34
    YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq    RSS    AIC
    <none>               33.370 19.341
    + EC     1    0.4063 32.963 20.949
    + Tair   1    0.1978 33.172 21.151
    + Tsoil  1    0.0801 33.289 21.264
    - MOIS   1   11.4629 44.832 26.790
    summary(model_Thur_Jun22)
    
    Call:
    lm(formula = YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M, 
        data = thur.jun.22)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -2.4725 -0.3853 -0.0512  0.5106  2.5927 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
    (Intercept)  1.25689    1.23111   1.021  0.31790   
    T           -0.70152    0.85285  -0.823  0.41920   
    S           -0.28573    0.85199  -0.335  0.74039   
    M           -1.33047    0.85187  -1.562  0.13198   
    MOIS         0.08097    0.02881   2.811  0.00992 **
    T:S          0.90545    1.21536   0.745  0.46382   
    T:M          3.16417    1.20452   2.627  0.01507 * 
    S:M          1.60516    1.20622   1.331  0.19632   
    T:S:M       -3.33147    1.70502  -1.954  0.06298 . 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 1.205 on 23 degrees of freedom
    Multiple R-squared:  0.4221,    Adjusted R-squared:  0.2211 
    F-statistic:   2.1 on 8 and 23 DF,  p-value: 0.07843

    Tillage, when combined with manure, showed a significant positive effect. Moisture seems to have a very high significance with a p-value of 0.00992 with a positive effect.

  4. Thurgem July 2023

    Variable DW was avoided because of perfect correlation with YIELD.

    im <- lm(YIELD ~ T * S * M, data = thur.july.23)
    
    model_Thur_Jul23 <- step(im,
                            scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + Tair + Tsoil + MOIS + Nconc + Cconc +
                            Block
                        ),
                        direction = "both")  
    Start:  AIC=120.51
    YIELD ~ T * S * M
    
            Df Sum of Sq    RSS    AIC
    + Cconc  1    375.54 463.02 103.50
    + Tair   1    107.83 730.74 118.11
    <none>               838.56 120.51
    + MOIS   1     14.73 823.83 121.94
    + Tsoil  1     11.83 826.73 122.06
    + Block  1      1.58 836.98 122.45
    + Nconc  1      0.16 838.40 122.50
    
    Step:  AIC=103.51
    YIELD ~ T + S + M + Cconc + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq    RSS    AIC
    <none>               463.02 103.50
    + Tair   1     25.27 437.75 103.71
    + Nconc  1     23.34 439.68 103.85
    + MOIS   1     13.93 449.09 104.53
    + Block  1      8.86 454.16 104.89
    + Tsoil  1      0.89 462.13 105.44
    - Cconc  1    375.54 838.56 120.51
    summary(model_Thur_Jul23)
    
    Call:
    lm(formula = YIELD ~ T + S + M + Cconc + T:S + T:M + S:M + T:S:M, 
        data = thur.july.23)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -6.3913 -3.3870  0.1161  3.1192  7.5971 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)  -17.675      6.975  -2.534 0.018544 *  
    T              9.336      3.417   2.732 0.011889 *  
    S              3.218      3.236   0.995 0.330326    
    M              6.331      3.577   1.770 0.089949 .  
    Cconc         21.942      5.080   4.319 0.000254 ***
    T:S           -9.041      4.551  -1.987 0.059001 .  
    T:M           -9.553      4.874  -1.960 0.062238 .  
    S:M           -5.955      4.826  -1.234 0.229726    
    T:S:M          5.002      6.663   0.751 0.460427    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 4.487 on 23 degrees of freedom
    Multiple R-squared:  0.5609,    Adjusted R-squared:  0.4082 
    F-statistic: 3.673 on 8 and 23 DF,  p-value: 0.006751

    The treatment tillage had a significant positive effect, whereas tillage combined with soil inoculum and manure had a significant negative effect. Manure also showed a significant positive effect, although its effect is not as pronounced as tillage, with an estimate of 6.331, whereas tillage had 9.336. Carbon concentration had a very significance with a p-value of 0.000254 positive effect.

  5. Thurgem Dec 2023

    im <- lm(YIELD ~ T * S * M , data = thur.dec.23)
    
    model_Thur_dec23 <- step(im,
                            scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + MOIS + P
                        ),
                        direction = "both")  
    Start:  AIC=10.88
    YIELD ~ T * S * M
    
           Df Sum of Sq    RSS    AIC
    <none>              27.267 10.878
    + MOIS  1   0.64431 26.622 12.113
    + P     1   0.46632 26.800 12.326
    summary(model_Thur_dec23)
    
    Call:
    lm(formula = YIELD ~ T * S * M, data = thur.dec.23)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -1.6688 -0.6007 -0.1421  0.6171  2.6133 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)   2.4113     0.5329   4.524 0.000139 ***
    T             1.0717     0.7537   1.422 0.167918    
    S             0.6397     0.7537   0.849 0.404415    
    M             0.8078     0.7537   1.072 0.294468    
    T:S          -0.5609     1.0659  -0.526 0.603575    
    T:M          -0.7722     1.0659  -0.724 0.475789    
    S:M           0.2662     1.0659   0.250 0.804918    
    T:S:M        -0.5409     1.5074  -0.359 0.722870    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 1.066 on 24 degrees of freedom
    Multiple R-squared:  0.1976,    Adjusted R-squared:  -0.03645 
    F-statistic: 0.8443 on 7 and 24 DF,  p-value: 0.5624

    Measurement done after a month (July 2023) showed no significant effects for any of the treatments, and the model obtained was fairly poor, with an R-squared of 19.76%.

  6. Kapsorok June 2022

    im <- lm(YIELD ~ T * S * M , data = kaps.jun.22)
    
    model_Kaps_jun22 <- step(im,
                            scope = list(
                          lower = ~ T * S * M, 
                          upper = ~ T * S * M + MOIS + Tsoil + Tair + EC + Block
                        ),
                        direction = "both")  
    Start:  AIC=51.3
    YIELD ~ T * S * M
    
            Df Sum of Sq    RSS    AIC
    + Block  1   20.3776 76.072 45.710
    <none>               96.449 51.305
    + EC     1    5.3695 91.080 51.472
    + MOIS   1    0.8644 95.585 53.017
    + Tair   1    0.5081 95.941 53.136
    + Tsoil  1    0.0257 96.423 53.296
    
    Step:  AIC=45.71
    YIELD ~ T + S + M + Block + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq    RSS    AIC
    + MOIS   1   13.2794 62.792 41.571
    <none>               76.072 45.710
    + Tair   1    0.1182 75.953 47.660
    + Tsoil  1    0.0763 75.995 47.678
    + EC     1    0.0501 76.022 47.689
    - Block  1   20.3776 96.449 51.305
    
    Step:  AIC=41.57
    YIELD ~ T + S + M + Block + MOIS + T:S + T:M + S:M + T:S:M
    
            Df Sum of Sq    RSS    AIC
    <none>               62.792 41.571
    + EC     1     0.585 62.207 43.271
    + Tair   1     0.336 62.456 43.399
    + Tsoil  1     0.242 62.551 43.448
    - MOIS   1    13.279 76.072 45.710
    - Block  1    32.793 95.585 53.017
    summary(model_Kaps_jun22)
    
    Call:
    lm(formula = YIELD ~ T + S + M + Block + MOIS + T:S + T:M + S:M + 
        T:S:M, data = kaps.jun.22)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -2.9130 -0.5417 -0.0770  0.4879  4.2881 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
    (Intercept)  2.88600    1.34779   2.141  0.04358 * 
    T            2.15809    1.19620   1.804  0.08492 . 
    S           -1.10692    1.19598  -0.926  0.36473   
    M           -0.95527    1.19678  -0.798  0.43329   
    Block       -1.35670    0.40025  -3.390  0.00264 **
    MOIS         0.12675    0.05876   2.157  0.04220 * 
    T:S          1.18085    1.69415   0.697  0.49309   
    T:M          0.06027    1.73371   0.035  0.97258   
    S:M          3.45978    1.71023   2.023  0.05539 . 
    T:S:M       -4.98098    2.49194  -1.999  0.05813 . 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 1.689 on 22 degrees of freedom
    Multiple R-squared:  0.594, Adjusted R-squared:  0.4279 
    F-statistic: 3.577 on 9 and 22 DF,  p-value: 0.007063

    Measurements showed a significant negative effect for tillage combined with soil inoculum and manure, tillage showed a positive effect, and soi inoculum combined with manure showed a positive effect. Moisture showed a positive effect with a p-value of 0.04220.

    Testing the significance of Block on dependent variable

    summary(aov(YIELD ~ Block, data = kab.nov.22))
                Df Sum Sq Mean Sq F value Pr(>F)
    Block        1    9.0   8.978   0.745  0.395
    Residuals   30  361.7  12.056               
    summary(aov(YIELD ~ Block, data = kab.dec.23))
                Df Sum Sq Mean Sq F value Pr(>F)
    Block        1    1.8   1.801   0.304  0.585
    Residuals   30  177.5   5.917               
    summary(aov(YIELD ~ Block, data = thur.jun.22))
                Df Sum Sq Mean Sq F value Pr(>F)
    Block        1   2.58   2.576   1.401  0.246
    Residuals   30  55.17   1.839               
    summary(aov(YIELD ~ Block, data = thur.july.23))
                Df Sum Sq Mean Sq F value Pr(>F)
    Block        1    1.6    1.58   0.045  0.833
    Residuals   30 1053.0   35.10               
    summary(aov(YIELD ~ Block, data = thur.dec.23))
                Df Sum Sq Mean Sq F value Pr(>F)
    Block        1   1.43   1.426   1.314  0.261
    Residuals   30  32.56   1.085               
    summary(aov(YIELD ~ Block, data = kaps.jun.22))
                Df Sum Sq Mean Sq F value Pr(>F)  
    Block        1  20.38  20.378   4.552 0.0412 *
    Residuals   30 134.29   4.476                 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Except for Kapsorok Block did not have a significant effect on YIELD. The effect of Block for Kapsorok is significant with a p-value of 0.0412.


    Creating table to summarize results for report.

    stargazer(model_Kaps_jun22,
              type = "text", 
              title = "Comparison of Regression Models across datssets with YIELD as dependent
              variable")

Secondary analysis

Considering the variable count the datasets on Kabianga November 2022 was only chosen for principal component analysis.

Kabianga November 2023

Performing PCA and selecting principal components

pca_result.kab <- kab.nov.22 |>
  select(-`Plot No.`, -M, -S, -T, - Treatment, -Code, -Block) |>
  PCA(scale.unit = TRUE)

fviz_eig(pca_result.kab, addlabels = TRUE, ylim = c(0, 50))

Corelation table and biplot

var_coord <- pca_result.kab$var$coord[, 1:2]

selected_vars <- var_coord[c("Tsoil", "Tair", "MOIS", "YIELD", "NH4", "NO3", "P", "pH", "EC", "MBC", "MBN", "Nconc", "Cconc"), ]

correlation_table <- data.frame(Variable = rownames(selected_vars), round(selected_vars, 2))

corelation_tibble <- as_tibble(correlation_table) |>
  rename(PC1 = Dim.1,
         PC2 = Dim.2)


print(corelation_tibble)
# A tibble: 13 × 3
   Variable   PC1   PC2
   <chr>    <dbl> <dbl>
 1 Tsoil    -0.61  0.52
 2 Tair     -0.69  0.5 
 3 MOIS      0.33  0.71
 4 YIELD     0.2   0.74
 5 NH4       0.67  0.18
 6 NO3       0.72 -0.27
 7 P        -0.49  0.22
 8 pH       -0.02 -0.06
 9 EC        0.76 -0.24
10 MBC       0.75  0.41
11 MBN       0.75  0.31
12 Nconc     0.88 -0.08
13 Cconc     0.8   0.19

K-means clustering based on principal components

Result may vary from report because set.seed(1234567) was added after report was written.

pc <- data.frame(pca_result.kab$ind$coord[, 1:2])
fviz_nbclust(pc, FUNcluster = kmeans, method = "wss")

set.seed(12334567)
k_means <- kmeans(pc, centers = 4)

fviz_pca_biplot(pca_result.kab,
                habillage = as.factor(k_means$cluster),
                addEllipses = TRUE
               )

fviz_pca_biplot(pca_result.kab,
                habillage = as.factor(kab.nov.22$Code)
               )

p_cat <- tibble(
  Code = kab.nov.22$Code,
  Cluster = as.factor(k_means$cluster)
)

ggplot(p_cat, aes(x = Code, fill=Cluster)) +
  geom_bar(stat = "count") +
  theme_minimal()

It seems as though manure treatment is associated with improved soil properties which are positively correlated with YIELD, except for EC and NO3, which seem to be poorly correlated. Tilage seems to have a very little to no impact on YIELD.

Just a plot to visualise variables and their correlation with PC1 and PC2

correlation_data <- data.frame(
  Variable = rownames(var_coord), 
  PC1 = pca_result.kab$var$coord[, 1], 
  PC2 = pca_result.kab$var$coord[, 2]
)



# Correlation plot
ggplot(correlation_data, aes(x = PC1, y = PC2, label = Variable)) +
  geom_point(size = 3, color = "blue") +  
  geom_text(vjust = 1.5, color = "black") +  
  theme_minimal() +
  labs(title = "Correlations of Variables with Principal Components",
       x = "PC1 (Principal Component 1)",
       y = "PC2 (Principal Component 2)")

Assesing the effect of solar radiation flux and precipitation

Creating tibbble for analysis

wi_kab_nov_22 <- read_csv("data/Weather data/wi_kab_dec_23.csv")
wi_kab_dec_23 <- read_csv("data/Weather data/wi_kab_nov_22.csv")


wi_thur_jun_22 <- read_csv("data/Weather data/wi_thur_jun_22.csv")
wi_thur_jul_23 <- read_csv("data/Weather data/wi_thur_jul_23.csv")
wi_thur_dec_23 <- read_csv("data/Weather data/wi_thur_dec_23.csv")

wi_kaps_jun_22 <- read_csv("data/Weather data/wi_kaps_jun_22.csv")


a <- kab.nov.22 |>
  select(YIELD, Block, Code, MOIS) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_kab_nov_22$Precipitation),
         Mean_srf = mean(wi_kab_nov_22$solar_radiation_flux),
         date = "nov_22")


b <- kab.dec.23 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_kab_dec_23$Precipitation),
         Mean_srf = mean(wi_kab_dec_23$solar_radiation_flux),
         date = "dec_23")

c <- thur.jun.22 |>
  select(YIELD, Block, Code, MOIS) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_thur_jun_22$Precipitation),
         Mean_srf = mean(wi_thur_jun_22$solar_radiation_flux),
         date = "jun_22")


d <- thur.july.23 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_thur_jul_23$Precipitation),
         Mean_srf = mean(wi_thur_jul_23$solar_radiation_flux),
         date = "jul_23")

e <- thur.dec.23 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_thur_dec_23$Precipitation),
         Mean_srf = mean(wi_thur_dec_23$solar_radiation_flux),
         date = "dec_23")



weather_info_kab <- bind_rows(a, b)

weather_info_thur <- bind_rows(c, d ,e)

weather_info_kaps <- kaps.jun.22 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_kaps_jun_22$Precipitation),
         Mean_srf = mean(wi_kaps_jun_22$solar_radiation_flux),
         date = "jun_22")

rm(a,b,c,d,e)
  1. Kabianga

    ggplot(weather_info_kab, aes(x = Mean_precip, y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) +
      geom_point() +
      geom_text_repel(aes(label = date)) +
      facet_wrap(~Code) +
      theme_minimal() +
      labs(title = "Yield Trends Over Time by Treatment For Precipitation", x = "Mean precipitation", y = "Yield")

    ggplot(weather_info_kab, aes(x = log(Mean_srf), y = mean_yield, group = Code, color = Code)) +
      geom_line(aes(linetype = Code)) +
      geom_point() +
      geom_text_repel(aes(label = date)) +
      facet_wrap(~Code) +
      theme_minimal() +
      labs(title = "Yield Trends Over Time by Treatment For Solar radiation Flux ", x = "Mean solar radiation flux (log scale)", y = "Yield")

    summary(lm(mean_yield ~ Mean_precip:Code, data = weather_info_kab))
    
    Call:
    lm(formula = mean_yield ~ Mean_precip:Code, data = weather_info_kab)
    
    Residuals:
         Min       1Q   Median       3Q      Max 
    -1.59945 -1.16581  0.05666  0.87641  2.26297 
    
    Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
    (Intercept)          35.7426     2.5642  13.939 2.31e-06 ***
    Mean_precip:CodeC    -5.3420     0.6415  -8.328 7.05e-05 ***
    Mean_precip:CodeM    -5.2672     0.6415  -8.211 7.71e-05 ***
    Mean_precip:CodeMS   -5.1600     0.6415  -8.044 8.80e-05 ***
    Mean_precip:CodeMST  -5.6693     0.6415  -8.838 4.80e-05 ***
    Mean_precip:CodeMT   -5.4248     0.6415  -8.457 6.38e-05 ***
    Mean_precip:CodeS    -4.5853     0.6415  -7.148 0.000186 ***
    Mean_precip:CodeST   -5.1048     0.6415  -7.958 9.43e-05 ***
    Mean_precip:CodeT    -5.4642     0.6415  -8.518 6.09e-05 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 1.769 on 7 degrees of freedom
    Multiple R-squared:  0.9279,    Adjusted R-squared:  0.8454 
    F-statistic: 11.25 on 8 and 7 DF,  p-value: 0.002275
    summary(lm(mean_yield ~ Mean_srf:Code, data = weather_info_kab))
    
    Call:
    lm(formula = mean_yield ~ Mean_srf:Code, data = weather_info_kab)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -2.0193 -0.7584 -0.0061  0.7237  1.8153 
    
    Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
    (Intercept)      -6.176e+01  7.236e+00  -8.535 6.02e-05 ***
    Mean_srf:CodeC    3.680e-06  3.604e-07  10.211 1.86e-05 ***
    Mean_srf:CodeM    3.711e-06  3.604e-07  10.298 1.76e-05 ***
    Mean_srf:CodeMS   3.715e-06  3.604e-07  10.311 1.75e-05 ***
    Mean_srf:CodeMST  3.597e-06  3.604e-07   9.982 2.17e-05 ***
    Mean_srf:CodeMT   3.641e-06  3.604e-07  10.103 2.00e-05 ***
    Mean_srf:CodeS    3.849e-06  3.604e-07  10.681 1.38e-05 ***
    Mean_srf:CodeST   3.711e-06  3.604e-07  10.299 1.76e-05 ***
    Mean_srf:CodeT    3.636e-06  3.604e-07  10.090 2.02e-05 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 1.538 on 7 degrees of freedom
    Multiple R-squared:  0.9454,    Adjusted R-squared:  0.8831 
    F-statistic: 15.16 on 8 and 7 DF,  p-value: 0.0008942
  2. Thurgum

    ggplot(weather_info_thur, aes(x = Mean_precip, y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) +
      geom_point() +
      geom_text_repel(aes(label = date)) +
      facet_wrap(~Code) +
      theme_minimal() +
      labs(title = "Yield Trends Over Time by Treatment For Precipitation", x = "Mean
           precipitation", y = "Yield")

    ggplot(weather_info_thur, aes(x = log(Mean_srf), y = mean_yield, group = Code, color = Code)) +
      geom_line(aes(linetype = Code)) +
      geom_point() +
      geom_text_repel(aes(label = date)) +
      facet_wrap(~Code) +
      theme_minimal() +
      labs(title = "Yield Trends Over Time by Treatment For Solar radiation Flux ", x = "Mean
           solar radiation flux (log scale)", y = "Yield")

    summary(lm(mean_yield ~ Mean_precip:Code, data = weather_info_thur))
    
    Call:
    lm(formula = mean_yield ~ Mean_precip:Code, data = weather_info_thur)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -4.2480 -2.9100 -0.0426  2.0472  6.8268 
    
    Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
    (Intercept)            9.370      1.516   6.180 1.76e-05 ***
    Mean_precip:CodeC     -2.885      1.716  -1.681    0.114    
    Mean_precip:CodeM     -2.837      1.716  -1.653    0.119    
    Mean_precip:CodeMS    -2.164      1.716  -1.261    0.227    
    Mean_precip:CodeMST   -3.177      1.716  -1.851    0.084 .  
    Mean_precip:CodeMT    -2.216      1.716  -1.291    0.216    
    Mean_precip:CodeS     -2.639      1.716  -1.537    0.145    
    Mean_precip:CodeST    -2.689      1.716  -1.567    0.138    
    Mean_precip:CodeT     -2.223      1.716  -1.295    0.215    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 3.887 on 15 degrees of freedom
    Multiple R-squared:  0.3227,    Adjusted R-squared:  -0.03853 
    F-statistic: 0.8933 on 8 and 15 DF,  p-value: 0.5452
    summary(lm(mean_yield ~ Mean_srf:Code, data = weather_info_thur))
    
    Call:
    lm(formula = mean_yield ~ Mean_srf:Code, data = weather_info_thur)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -4.5003 -2.4308  0.3415  1.9385  4.8084 
    
    Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
    (Intercept)       7.197e+01  1.566e+01   4.595 0.000350 ***
    Mean_srf:CodeC   -3.230e-06  7.686e-07  -4.203 0.000769 ***
    Mean_srf:CodeM   -3.251e-06  7.686e-07  -4.229 0.000728 ***
    Mean_srf:CodeMS  -3.182e-06  7.686e-07  -4.140 0.000873 ***
    Mean_srf:CodeMST -3.310e-06  7.686e-07  -4.307 0.000623 ***
    Mean_srf:CodeMT  -3.169e-06  7.686e-07  -4.123 0.000903 ***
    Mean_srf:CodeS   -3.215e-06  7.686e-07  -4.183 0.000799 ***
    Mean_srf:CodeST  -3.242e-06  7.686e-07  -4.218 0.000745 ***
    Mean_srf:CodeT   -3.165e-06  7.686e-07  -4.118 0.000913 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    Residual standard error: 3.093 on 15 degrees of freedom
    Multiple R-squared:  0.571, Adjusted R-squared:  0.3422 
    F-statistic: 2.495 on 8 and 15 DF,  p-value: 0.06039

In Kabianga, productivity was lower during the months with high rainfall. However, when considering solar radiation flux, yield increased with higher levels of solar radiation. In contrast, at Thurgem, YIELD decreased as both rainfall and solar radiation flux increased.

Miscellaneous